Anomalous diffusion of driven particles in supercooled liquids 
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We have performed non-equilibrium dynamics simulations of a binary Lennard- Jones mixture in 
which an external force is applied on a single tagged particle. For the diffusive properties of this 
particle parallel to the force superdiffusive behavior at intermediate times as well as giant long-time 
diffusivity is observed. A quantitative description of this non-trivial behavior is given by a continuous 
time random walk analysis of the system in configuration space. We further demonstrate, that the 
same physical properties which are responsible for the superdiffusivity in non-equilibrium systems 
also determine the non-Gaussian parameter in equilibrium systems. 



Introduction. Due to the distinct multi-particle dy- 
namics of glass-forming systems several interesting prop- 
erties can be observed like the occurrence of dynamical 
heterogeneities [T] or the violation of the Stokes-Einstcin 
relation [U [3] . In the non-equilibrium situation, the ob- 
served phenomena can become even more complex. Re- 
cently, Winter et al. performed computer simulations of 
a tracer particle which is driven by a constant external 
field trough a binary Yukawa fluid [4] . It was shown that 
for this microrheological simulation the diffusive proper- 
ties of the particle become highly anisotropic: While the 
mean squared displacement (MSD) of the tracer parti- 
cle perpendicular to the force direction (a;^)(f) increases 
with increasing force but still displays a diffusive behav- 
ior, the centered MSD parallel to the force direction 

a\t) = {xl{t))-{x\\{t)f (1) 

displayed a superdiffusive behavior at the observed 
time range. This result was rationalized in terms of a 
special type of biased trap model [5] in which a superdif- 
fusive behavior is predicted due to rising fluctuations. 
Therefore, it was stated that the diffusion constant for 
the parallel direction of the tracer particle does not exist 
[3]. However, this model has to be regarded with care be- 
cause the rising fluctuations would lead to a permanently 
increasing energy barriers. This scenario is difficult to 
reconcile with the observed stationary behavior. Mode- 
coupling theory has been very successful to predict, e.g., 
the nonlinear mobility dependence on the applied force 
in microrheological simulations [5HSj- Interestingly, this 
superdiffusive behavior could not be reproduced [S]. 

A different approached is used by Jack et al. [TU] . Mo- 
tivated by the analysis of an one-dimensional spin facil- 
itated models, they performed an analytical calculation 
for a biased continuous time random walk (CTRW). For 
this model, a diffusive regime is predicted for long times. 
This diffusive regime is characterized by a strong depen- 
dence on the width of the used waiting time distribution. 
Broader waiting time distributions lead to a dramatic 
increase of spatial fiuctuations, denoted as "giant diffu- 
sivity" [TT]- 



The key goal of this paper is to elucidate the prop- 
erties of the superdiffusivity in the driven particle dy- 
namics. First, we present a formal expression which re- 
lates the superdiffusivity to dynamic heterogeneities in 
the CTRW framework. Second, for the trajectories of a 
glass-forming model system we can extract the relevant 
observables from an appropriate CTRW analysis and pre- 
dict the superdiffusive behavior in a quantitative way. 
For the long-time limit our expression reduces to the gi- 
ant diffusivity as calculated in Ref.[Tn]. Third, we are 
able to show that the superdiffusivity has a deep physi- 
cal connection to the non-Gaussian parameter (NGP) in 
equilibrium, thereby establishing a strong connection be- 
tween the non-equilibrium and the equilibrium dynamics 
of glass-forming systems. 

Simulations. We have performed computer simula- 
tions of a binary mixture of Lennard- Jones particles [2] 
(BMLJ) which we have extended by applying a constant 
force on one randomly selected particle. Constant tem- 
perature conditions are ensured by using a Nose-Hoover 
thermostat |13j . By applying a suitable minimization 
procedure it is possible to track minima of the PEL, 
called inherent structures (IS), which the system had ex- 
plored during its time evolution. Analogous to equilib- 
rium simulations |14H16| we have recently demonstrated 
|17j . that the time evolution of a small stationary non- 
equilibrium system (consisting of 65 particles, denoted 
as BMLJ65) can be analogous to equilibrium systems 
[m [15] described in terms of a continuous time random 
walk (CTRW) of the system between coarse grained min- 
ima, called metabasins (MB). This projection allows for 
a discretization of the system trajectory into dynami- 
cal events, which are characterized by the distribution 
of particle displacements during one transition, and a 
broad waiting time distribution. As shown in |17j the 
linear and nonlinear response only shows very small fi- 
nite size effects. Here we will also show that the results 
of this work can be transferred to the properties of large 
systems as well. 

Results. Focusing on the diffusive behavior of the 
tracer particle parallel to the force direction, our ap- 
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FIG. 1. Centered MSD cr^(n) of BMLJ65 as a function of the 
number of transitions n at a temperature T — 0.475. The 
dashed hues indicate the diffusive lengths a| parallel to the 
force direction. 



proach offers two different routes to define the centered 
MSD: On the one hand, one can consider the centered 
MSD after a certain number of MB transitions n on the 
other hand, it can be evaluated after a certain time t. In 
the following we will distinguish between these quanti- 
ties by writing (T^(n) and cr^(t), respectively. Similar to 
the equilibrium dynamics |14j , (n) grows linearly after 
more than ^ 20 transitions (see Fig. [ij . In marked con- 
trast, cr^(i) displays a superdiffusive behavior (see Fig.[2| 
as it was reported for the binary Yukawa fluid 4 • From 
(n) one can define for large n the diffusive length scale 
a?, via 



= lim 



(2) 



see also PTI. 

To achieve a more quantitative understanding and 
to unravel the surprising qualitative differences between 
cr^(ri) and cr^(t), we have performed an analytical calcu- 
lation of cr^(i) within the CTRW framework. We start 
with a one-dimensional CTRW with an elementary step 
^i.W — + ^^\\- ^^\\ denotes the average translation 
the particle performs along the force direction during one 
MB transition. As shown in |I7j, Axy is basically propor- 
tional to F in the whole force interval considered in this 
work, flj II is considered to be the remaining translational 
length with (an) = 0. Successive steps are regarded as 



uncorrelated so that 



For reasons of 



consistency with our previous work, we will further de- 
note {ap simply as Oy. Then, the MSD of the particle is 
given by the sum over all steps n which were performed 
up to a time t: 



n(t) 

(4w)-((E(«-ii+^^ii))' 



(3) 



This yields for the time evolution of the MSD 

{xl{t))^{n{t))a\ + {n\t))Ax\. (4) 



FIG. 2. Centered MSD a'^{t) of BMLJ1560 as a function of 
time t at a temperature T = 0.475. The dashed lines corre- 
spond to the theoretical prediction by Eq. [5] The gray bars 
at short times indicate the equilibrium-like diffusion (Eq. [6]l , 
at long times the long-time diffusion (Eq. |I0[ ) under consider- 
ation of the numerical error. 



In Eq. |4] (n) denotes the average number of jumps the 
particle performs in a certain time t and (n^) the fluctu- 
ation of these, respectively. Considering cr^(t) instead of 
the MSD one has to subtract the squared first moment 
of the particle displacement which is given by(x|[(t))^ = 

(n(i))'Axf. 

After subtracting this expression from Eq. |4]one finally 
obtains 



a2(t) = {nit)) a\ + [{n^{t)) - {n{t)f] A.Tjj 



(5) 



The first term of Eq. [5] can be identified with the one- 
dimensional equilibrium-like diffusive process 



2 Dh t = at 



(6) 



We use the term "equilibrium-like" because, as we have 
shown for our numerical data in |17j , both and (r) dis- 
play a force dependence so that the actual value of I?|| 
increases with increasing force. The superdiffusivity of 
CT^(t) can be related to the latter term which is only vis- 
ible in driven systems. The expression [(n^(t)) — {n{t)) ] 
describes the heterogeneity of the number of occurring 
jumps in a certain time interval and can be directly ob- 
tained from our trajectories. 

For BMLJ65 one is able to perform a complete CTRW 
analysis so that each observable in Eq. [5] is directly ac- 
cessible. As it can be seen in Fig. [2] this ansatz allows to 
quantitatively reproduce the behavior of cr^(t). 

The time evolution can be divided into three regimes: 
At short times one observes a decay of cr^(i) due to lo- 
cal caging until a short diffusive regime is reached at 
t « (r) which reflects the equilibrium-like diffusion pro- 
cess (Eq. |6| . The initial decay of cr^ (t) /t can be qualita- 
tively understood by considering, that a'^{n)/n displays 
a similar decay due to slightly forward-backward corre- 
lations until the constant diffusive length of, is reached 
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at n w 20 transitions (see [2] for further details). It is 
important to notice that only in case of small forces the 
minimum of a'^{t)/t indicates the true value of Dy while 
at higher forces it is already superimposed by superdiffu- 
sive contributions. At intermediate times, one observes a 
superdiffusive behavior which is caused by the nonlinear 
evolution of [{■n?{t)) — {nit))'^]. At long times, when t is 
larger than the largest measured waiting time, one indeed 
observes an indication, that the MSD becomes diffusive 
again but with a significantly larger diffusion constant. 
For this particular long time behavior one is able to yield 
an analytical prediction by the present CTRW ansatz: 

The waiting time distribution tf{T) of the CTRW can 
be characterized by its average value (r) and its variance 
V = (r^) — (r^). Due to the central limit theorem the 
distribution of the cumulated waiting time t„ of a large 
number of jumps n, P„(r„), is given by 



lim P„(r„) cx exp 



{rn-n{T)f 
2Vn 



(7) 



The corresponding probability to find n jumps in a 
large time interval t, Pt{n), is directly related to P„(t„). 
With the substitution n ~ t/{T) and identifying r„ = t 
we obtain from Eq. [7] the expression 



lim Pt(n) cx exp 



'in- 



t 

t)) 



(8) 



Determination of the second moment of Pt{n) yields 



V 

t (r) 



- 1 



(9) 



and hence, by combining Eq. [9] and Eq. [5] for the long- 
time behavior of cr^(t) 



lim — = Duf = D« 

t^oo t " " 



- 1 



(10) 



In this equation, / describes the factor which re- 
lates the equilibrium-like and the long-time diffusion con- 
stants. Independent from our derivation Jack et al. ana- 
lytically obtained a similar result for the giant diffusivity 
by considering the Montroll- Weiss equation of a biased 
CTRW [To]. The long-time diffusion constant Dy/ was 
already indicated in Fig. [2] For BMLJ65, it is possible 
to explicitly compute the long-time diffusivity because 
of the direct access to the CTRW observables in Eq. [TO] 
Importantly, it is also possible to estimate the long-time 
behavior for larger system sizes. As shown in the sup- 
plementary material the numerically observed degree of 
superdiffusivity is fully compatible with the theoretical 
expectation. 

The presented ansatz allows one to give an explicit cri- 
terion how long a particle requires to reach the diffusive 
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FIG. 3. One-dimensional NGP 02 (t) multiplied by time t as a 
function of t. The dashed lines corresponds to the theoretical 
prediction in Eq. |12[ the arrows indicate the structural relax- 
ation time Ta for the different temperatures. Inset: 02 (t) t 
at a temperature T — 0.475 together with its temporal and 
spatial contributions (see text). 



regime. It is related to the applicability of the central 
limit theorem and thus to the width of the waiting time 
distribution: The narrower the waiting time distribution, 
the earlier the particle becomes diffusive. Since the ap- 
plication of a strong microrhological perturbation causes 
a narrowing of the waiting time distribution |17| one ex- 
pects an earlier advent of the long-time diffusivity at high 
forces. This behavior can be qualitatively observed in 
Fig.|2]as weh. 

Besides the MSD of a driven particle, the heterogeneity 
of MB transitions can also be observed in equilibrium sys- 
tems by analyzing the NGP a2 (t) of the one-dimensional 
particle displacement which is defined as 



a2(t) 



3{x^t)f 



(11) 



Using the ansatz of an unbiased CTRW one obtains 
for the NGP 

c..(0- [^^^^^f',y^^% /^nP.(n)A(.). (12) 
{n{t)) J 

The latter term describes the non-Gaussianity of the 
cumulated displacements a„ after n transitions 



Ain) 



3(a2) 



2 



(13) 



weighted by the probability to find exactly n transitions 
at a time t. In what follows we use the approximation 
that / dn Pt{n)A{n) « A{{n(t))). 

Eq. [12] contains two contributions: The first term in- 
cludes the heterogeneity of the performed jumps n in a 
certain time interval. It is the same quantity which was 
observed to be responsible for the superdiffusive behav- 
ior in the non-equilibrium system. Because this term is 
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independent of any length scales, one can regard it as 
a measure for the temporal heterogeneities of the sys- 
tem dynamics. The latter term reflects spatial hetero- 
geneities of the elementary MB transition which become 
less important at a larger number of transitions because 
the distribution of the cumulated lengths approaches a 
Gaussian shape. 

In Fig.|3]a2(i) -t is shown at different temperatures to- 
gether with the theoretical prediction by Eq. [T2j Please 
note that, as it was also shown by Liao et al. [TB], in 
case of transitions between adjacent MB, a2{t) displays 
a monotonic decay. This behavior can be understood 
by considering that the initial growth of a2 (t) is caused 
by vibrational parts (short times) and the /3-relaxation 
process (at intermediate times) [HJ [inj while, by con- 
struction, MB trajectories only resolve the a-relaxation 
process [2D|. As one can see, the theoretical predic- 
tion allows one to fully describe the NGP at each tem- 
perature. One further observes that for long times, 
a2{t) ■ t approaches a constant which corresponds to a 
decay of a2 (i) oc 1/t which is exactly the expectation for 
[(7i^(t)) — {n{tyf']/ {n{t))'^ when the central limit theorem 
becomes valid (see Eq. [o]) . 

In the inset of Fig. [3] the different temporal and spatial 
contributions to a2{t) ■ t are shown. At very short times, 
the behavior of a2{t)-t is mainly determined by A{{n{t))) 
while above t « 10^, the temporal part of Eq. [l2]is found 
to be the major contribution. At t « 10^, one observes 
A{{n{t))) a 1/t which indicates that the central limit 
theorem starts to hold for the distribution of the spatial 
displacement. Indeed, for the waiting time distribution 
the central limit theorem is only fulfilled on a larger time 
scale, so that there is still a growth of a2{t) ■ t. 

It is know from a comparison between mode-coupling 
theory and Brownian dynamics simulations of BML J [5T] , 
that mode-coupling theory tends to strongly underesti- 
mate the magnitude of a2{t) in the diffusive regime. This 
differences between theory and simulation are also known 
for the hard-sphere system [2^ Recently, it was 

further reported, that simplified mode-coupling theory 
models cannot reproduce the superdiffusive behavior of 
a driven particle along its force direction [9]. Therefore, 
it is quite remarkable that the CTRW approach enables 
to relate both, the non-Gaussianity of the equilibrium 
system and the superdiffusive behavior of the stationary 
non-equilibrium system to have the same origin, reflect- 
ing the presence of the dynamic heterogeneities. One 
thus might argue whether mode-coupling theory, possi- 
bly due to its dependence on average quantities |22| , is 
not able to fully describe these fluctuations. As a conse- 
quence both effects cannot be qualitatively reproduced. 

Summary. In the present paper we have demonstrated, 
that a model of a biased CTRW allows to fully predict 
the anomalous diffusion of a driven particle in a super- 
cooled medium which is characterized by equilibrium- 
like diffusion, superdiffusivity and long-time diffusivity. 



It was further shown, that the origin of the superdiffu- 
sivity results from temporal fluctuations of the system 
dynamics which become visible due the applied bias. 
Indeed, also in equilibrium the same fluctuations are 
present and determine the evolution of the NGP a2{t)- 
Therefore, the connection between superdiffusive behav- 
ior and non-Gaussianity is a remarkable example, how 
non-equilibrium dynamics also enables a deeper physical 
understanding of the equilibrium system, by uncovering 
essential underlying physical properties. 
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more, C. F. E. Schroer thanks the NRW Graduate School 
of Chemistry for funding. We acknowledge helpful dis- 
cussions with J. Horbach, C. Rehwald and D. Winter. 
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LONG-TIME DIFFUSIVITY OF LARGER 
SYSTEMS 

As we had demonstrated, the continuous time random 
walk (CTRW) analysis allows one to predict the long- 
time diffusivity of the driven particle as 



lim 



t 



= D« f = D« 




- 1 



(14) 

For a binary mixture of = 65 Lennard- Jones parti- 
cles (BMLJ65) the required CTRW quantities can be ex- 
plicitly computed from the metabasin (MB) trajectories. 
For systems larger then BMLJ65 it is not possible be- 
cause there is no direct access to the required quantities. 
Unfortunately it is known, that especially the structural 
relaxation time r^, which reflects the higher moments of 
the waiting time distribution, displays significant finite- 
size effects [ini m]. Therefore, one must expect, that 
the degree of superdiffusivity will be different at larger 
system sizes. It is nevertheless possible to estimate the 
magnitude of long-time diffusivity from the real space 
trajectories of larger systems. 

We start with a relation between the average waiting 
time (t) and the diffusion constants ^||/_l of the tracer 
particle parallel and perpendicular to the direction of the 
force. It is given by 



II/-L 



(15) 



,j_ denote the apparent diffusive lengths dur- 



where 

ing one MB transition [HIITT]. Furthermore, the second 
moment of the waiting time distribution (r^) is connected 
to r„ by [UlEl] 



lim Ta {q) = Ta 

q—>OQ 



2(r) 



(16) 



Please note, that in a strict sense the relation between 
liniq^oo Ta{q) and {t"^) is only valid for BMLJ65. We will 
identify the relevant value of Tq, further below. Inserting 
Eq. 



15 and Eq. 16 in Eq. 14 yields 



a-, 



4 T„ 



- 1 



(17) 



In contrast to Eq. 14 Eq. 17 only depends on micro- 
scopic lengths and quantities which are measurable in 
real space. Please note, that the diffusion constant D|| 

was substituted by -^D±^ because is not accessible 
due to the superdiffusive behavior of cr(t). 
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FIG. 4. Ratios of the diffusion coefficient _Dx perpendicular to 
the force direction and the drift velocity v of BMLJ1560 and 
the corresponding microscopic length scales a\ and Aa;|| of 
BMLJ65 versus the applied Force F. Both systems were sim- 
ulated at T = 0.475. 



We now assume that the microscopic length scales in 
Eq.[T7]are the same for BMLJ65 and larger systems. This 
assumption is made because the microscopic lengths refer 
to the local spatial displacement of single particles during 
single relaxation processes which should be largely inde- 
pendent from the system size. To verify this, we have 
computed the ratio of the drift velocity v and the perpen- 
dicular diffusion constant D± of a system with A^ — 1560 
particles (BMLJ156G). Because v = ^ and Eq.[l5](for 
BMLJ65), this ratio can be identified with a^/2Ax|j. If 
the assumption holds, a^/Axy of BMLJ65 should be es- 
sentially the same as 2D±/v of BMLJ1560. The corre- 
sponding plot is shown in Fig. [4] One can observe, that 
the mismatch between these curve is less than 20% so 
that the independence of the microscopic length scales 
on system size indeed seems to be reasonable. 



The missing observable Tq, in Eq. 17 is determined by 
computing the incoherent scattering function 



S^{q,t) = (cos {q{x^{t)~x^{to)f)) 



(18) 



of the tracer particle perpendicular to the force direction. 
The reason for the choice of the perpendicular direction 
is, that one would expect, analogous to the superdiffu- 
sive behavior of cr^(i), an additional relaxation along the 
parallel direction which does not reflect the underlying 
waiting time distribution. Following the procedure of 
Rehwald et al. [21] we have computed the apparent wave- 
vector relaxation 



r{q) 



poo 

/ dtS±{q,t). 



(19) 



and fitted the resulting curves in the range 9 < 8 by 

1 



q^ ■ Da 



(20) 
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FIG. 5. Centered MSD a'^{t) versus time t for BMLJ1560 at 
T — 0.475. The solid lines at short times correspond to the 
equilibrium-like diffusion constant (Eq. 151, the solid lines at 
long times to the predicted long-time diffusivity (Eq |17[ ). 



finally allowing us to determine the local relaxation time 



Ta- 



in Fig.[5]cr2(t) of BMLJ1560 is shown together with the 
predicted equilibrium-like and long-time diffusion con- 
stants. One observe a quantitative agreement between 
the numerical data and the theoretical prediction. Com- 
pared to BMLJ65, the factor between equilibrium-like 
and long-time diffusion is smaller and the long-time dif- 
fusivity appears earlier. This behavior indicates a nar- 
rowing of the underlying waiting time distribution of 
BMLJ1560. One can rationalize this by regarding the 
large system as a composition of small elementary system 
which are coupled so that relaxation processes of single 
subsystems induce further relaxation events in adjacent 
systems. This concept has been discussed in ^(71 124j . 
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